Mott Physics and Topological Phase Transition in Correlated Dirac Fermions 
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We investigate the interplay between the strong correlation and the spin-orbital coupling in the 
Kane-Mele-Hubbard model and obtain the qualitative phase diagram via the variational cluster 
approach. We identify, through an increase of the Hubbard U, the transition from the topological 
band insulator to either the spin liquid phase or the easy-plane antiferromagnetic insulating phase, 
depending on the strength of the spin-orbit coupling. A nontrivial evolution of the bulk bands in 
the topological quantum phase transition is also demonstrated. 
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In recent years, a new field has emerged in condensed 
matter physics, based on the realization that a spin-orbit 
interaction can lead to topologically insulating electronic 
phases [HE]- A topological band insulator (TBI) has 
a nontrivial band structure resulting from the strong 
spin-orbit coupling. Theoretical and experimental stud- 
ies have found such materials in both two (2D) [3H5] and 
three (3D) [rJHTU] dimensions. A common property of 
TBI is that it has a charge excitation gap in the bulk, 
but with gapless helical edge (or surface) states protected 
by the time reversal symmetry lying inside the bulk in- 
sulating gap. As a new quantum state, which is the Zi- 
graded topological distinction from other conventional in- 
sulators, it has attracted great attention. Though great 
progress has been achieved, the current researches mostly 
focus on the weakly interacting systems. It has been pro- 
posed that the topological insulator may also appear in 
the systems with substantial electron correlations, such 
as Ad and 5d transition metal oxides [TT] [T2J • And the 
electron interaction effect plays a crucial role in determin- 
ing the ground state of topological insulators in the 2D 
limit [13] ■ Therefore, the effects of electron correlations 
on the topological insulators present a new challenge. 

The correlation effects in topological insulators can be 
studied either by interaction-driven topological insula- 
tors [MHIB] or by introducing interactions to a system 
with a strong spin-orbit coupling [ITJ [TT1 [18]. In this 
Letter, we investigate the model proposed by Kane and 
Mele [3] on the honeycomb lattice for describing a 2D 
topological insulator, and introduce the Hubbard inter- 
action to this model to analyze the Mott physics. Re- 
cently, the Hubbard model on the honeycomb lattice 
have been studied by Meng et al [19] using the quan- 
tum Monte Carlo (QMC) method, in which a spin liquid 
(SL) phase is found to exist between the semi-metallic 
(SM) phase and the antiferromagnetically (AF) ordered 
Mott insulator (MI) phase for a range of the on-site inter- 
action U . The mean field analysis and QMC simulations 
for the Kane-Mele-Hubbard (KMH) model reveal that 



the TBI phase is unstable against the magnetic order- 
ing phase [T71 [201 HI] • But the whole phase diagram of 
the KMH model, especially the transition between the 
TBI and the MI, and the nature of the single-particle 
excitations in the bulk and on the edges are still open 
theoretical questions. As the existence of gapless edge 
states is the direct manifestation of the topological na- 
ture, the study of the single-particle excitation spectra is 
the natural way to investigate the phase transition be- 
tween TBI and MI. Here, we use the (zero temperature) 
variational cluster approach (VCA) [22], which goes be- 
yond the mean field theory and takes into account ex- 
actly the effects of short-range correlations by an exact 
diagonlization of the separative clusters. We find a topo- 
logical quantum phase transition from TBI to MI with 
increasing U and this process shows a nontrivial evolu- 
tion. Starting from TBI, the spin-orbit coupling gap Aso 
closes first and then the Mott gap opens up but without 
the gapless edge states for increasing U, which is closely 
related to the topological properties of the system. The 
closing process of A50 driven by the correlations is ac- 
companying with a splitting of both the conduction and 
valence bands. In the strong spin-orbit coupling regime, 
the state transiting from TBI is the easy-plane AF Mott 
insulator. In the weak coupling regime, a spin liquid 
phase emerges between the TBI and the AF Mott insu- 
lators. In addition, we also find a decrease in the velocity 
of the helical edge states due to the correlations in the 
TBI phase. 

The Kane-Mele-Hubbard model is defined as H — Ho+ 
Hjj, where Ho is the model proposed by Kane and Mele 
on the honeycomb lattice as shown in Fig. [ija) [3J, 
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and Hjj the Hubbard interaction, 

H v = 11^11^. (2) 
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FIG. 1: (color online) (a) The 6-site cluster tiling (dashed 
lines) on honeycomb lattice used in our calculations of the 
bulk properties with VCA. A and B denote the two inequiv- 
alent sites, ai and a2 are the lattice unit vectors, (b) The 
first Brillouin zone, bi and b2 are the reciprocal-lattice vec- 
tors, (c) Schematic view of tiling the ribbon for calculating 
the edge states. Here, two superlattices (rectangle with solid 
lines ) are shown. Each superlattice makes up of two clusters 
containing 12 sites as separated by the dashed line). 

Here, and ((ij)) denote the nearest neighbor (NN) 
and the next nearest neighbor (NNN), respectively. A is 
the spin-orbit coupling constant and r the Pauli matrices. 
Vij = +1( — 1) if the electron makes a left (right) turn to 
get to the NNN site. Others are in standard notation. 

The VCA is a cluster method of the self-energy func- 
tional approach (SFA) [55], which approximates the self- 
energy of the original system by the self-energy X' of 
an exactly solvable reference system with the same in- 
teraction term. It has been successfully applied to, for 
instance, the problem of competing phases in high-T c su- 
perconductors [531 [53] . Despite the considerable finite- 
size errors, the VCA can predict the qualitatively correct 
trend for the phase diagram [25 . In VCA, the lattice is 
tiled into identical clusters (as illustrated in Fig. [I] each 
cluster contains a hexagon in our calculations) and the 
reference system is made up of the decoupled clusters. 
The single-particle parameters (denoted by t') of the ref- 
erence system are optimized according to the variational 
principle. And one can add any Weiss field to study the 
symmetry broken phases. For any X' parameterized as 
X'(t'), we have the grand potential: 

n[S'(t')] = ft'(t') +Trln[-(G - 1 - S'(t'))" 1 ] 

- Trln[-G'(t')], (3) 

where O'(t') and G'(t') are the grand potential and 
Green's function of the reference system, Go is the 
free Green's function without interactions. The phys- 
ical self-energy S is given by the stationary point 
dn[E'(t')]/dt' = 0. For any t', S'(t') is related to 
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FIG. 2: (color online) Qualitative phase diagram of KMH 
model. SM, TBI, SL and AF insulator denote the semi-metal, 
topological band insulator, spin liquid and antiferromagnetic 
insulator, respectively. Above the white dashed line in the 
AF insulator phase, the z-term of the AF order disappears. 

the lattice Green's function G(t') by the Dyson equation 
G _1 (t') = Gq 1 - S'(t'). G(t') can be determined via 
the cluster perturbation theory [26 , in which G'(t') is 
calculated by the exact diagonalization method and the 
intercluster hopping V is treated perturbatively. In mo- 
mentum space, G(t') can be expressed in terms of G'(t') 
and V as G(k,w) = G'(k, u)[l - V(k)G'(k, w)]" 1 . 

To calculate the edge states, we consider a strip geom- 
etry and construct a supercluster which is made of sev- 
eral clusters. In Fig. [ljc), for example, we arrange two 
clusters (12 sites in one cluster) in y direction to form a 
supercluster. The Green funciton of the supercluster is 
given as (G sc ) _1 = (G') _1 — W and W is the intercluster 
hopping matrix in the supercluster. 

To test the existence of the possible AF order, we will 
include the following Weiss field, 

Hi F = h a AF j2(- i r4Aci', (4) 

i 

where rji = or 1, when i £ A or B. In the absence of the 
spin-orbit interaction, the spin sector has a SU(2) sym- 
metry. So, we have h z AF = h x AF . However, this relation 
is broken down when the spin-orbit interaction is turned 
on. In this case, we will calculate the grand potential 
Q(/iaf) as a function of h AF and h AF , respectively. 

Our main results on the interplay between the Hub- 
bard interaction and the spin-orbit coupling are summa- 
rized in the U — A phase diagram [Fig. [2] . Let us first 
discuss the A = line. In VCA, the existence of the AF 
order can be determined by the h AF dependence of the 
grand potential Q(h A p). Fig. [3]Ja) presents the results 
for different Hubbard interactions U. For weak U, such 
as U = 2t and At, fl(h A p) shows a monotonic increase 
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FIG. 3: (color online) (a) f2 as a function of Haf for various 
values of U at A = 0. (b) The density of states for U — 4 and 
A = 0. (c) and (d): vs Yiaf at A = 0.2t along the z and 
^-directions, respectively. 



with h z AF (h AF — h x AF in this case), indicating that no 
AF order forms in the system. However, for a large U 
such as U > 64, a minimum appears at finite h z AF and 
this minimum moves to lower h z AF values with increase 
of U. Therefore, we can infer that an AF order exists for 
a large U as expected. Interestingly, when plotting the 
density of states (DOS) for U = At as shown in Fig. [3jb) , 
we find that an obvious Mott gap has opened up around 
the Fermi energy. This paramagnetic insulating phase is 
identified as the SL phase as also been found recently by 
Meng et al using the QMC simulation [T5]. Therefore, 
the system will undergo phase transitions from the semi- 
metal(SM) to SL and then to AF Mott insulator with 
the increase of U. Thus, we can reproduce the QMC 
simulation results calculated for A = |19j . 

When turning on the spin-orbit coupling, we find that 
the SL phase maintains for a range of spin-orbit coupling 
up to A = 0.1254. On the other hand, the AF order is not 
isotropic. As seen from Figs. |3jc) and (d), no minimum 
is found at U — At for A = 0.24 in the h z AF dependence, 
but it can be found in the h x AF dependence. It indi- 
cates that within a range of U, the z-direction AF order 
is destroyed once the spin-orbit coupling is present. For 
A < 0.254, when increasing U further, we find the appear- 
ance of the z-term in the AF order eventually. However, 
for A > 0.254, it has not been found up to U — 104. Thus, 
in the phase diagram we plot the white dashed-line sep- 
arating the AF order with and without the z-term. The 
easy-plane AF order is the result of the interplay between 
the Hubbard interaction and the spin-orbit coupling. As 
is well known, the NN hopping will generate an isotropic 
AF Heisenberg term Hi = JiJ2(ij) Sj • Sj with J\ = 
At 2 /U in the strong-coupling limit. Similarly, the NNN 
spin-orbit coupling generates an anisotropic exchanging 




(9) 








1 






1 K 1 


* r 




term H 2 = J 2 Y, m) (SfS : 



S?S] 



S Z S Z ) HE, with 



■i 
■I 



FIG. 4: (color online) A(k, w) for single-particle excitations 
in the bulk [Figs. (a), (c), (e) and (g)] and in a ribbon with 
the zigzag edges [Figs.(b), (d), (f) and (h)] as illustrated in 
FigHTc) at A = 0.14. The white dashed curves in Figs.(c), (e) 
and (g) are the mean-field fits discussed in the text. From 
the up to down figures, U — 0, 24, 3t, At. The inset shows the 
^/-dependence of the renormalized velocity of edge states at 
A = 0.14 and 0.24. 



J2 = A\ 2 /U. Notice that the z term in H 2 favors an- 
tiparallel alignment of the spins on the NNN sites, thus 
it will introduce a frustration to the NN AF correlation 
expressed by H±. On the other hand, the xy term in 
H 2 favors a ferromagnetic alignment, so no frustration is 
introduced. As a result, the H 2 term coming from the 
spin-orbit coupling will suppress the z-term of the AF 
order. 

At another limit U = 0, a TBI is expected to occur 
once the spin-orbit coupling is turned on [3]. The TBI 
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is characterized by gapless edge states protected by the 
bulk gap opened by the spin-orbit coupling in the single- 
particle spectrum. The spectral function of single parti- 
cles is given by A(k,w) = — 2ImG(k,a;)/7r. The results 
for several U at A = O.lt are presented in Fig. |4j where 
the bulk bands are plotted along the lines shown in Fig. 
[jjb) and the edge states are calculated from a ribbon with 
the zigzag edges [Fig. [ijc)]. For [7 = 0, one can see that 
a bulk gap opens resulting from the spin-orbit coupling 
[Fig. |4ja)]. At the meantime, clear gapless edge states 
with sizeable spectral weights emerge [Fig. gb)]. These 
results reproduce perfectly the characters of a TBI [3]. 
In the presence of Hubbard interaction, the VCA cap- 
tures exactly the short-range correlation effects by the 
exact diagonlization of the small clusters used to tile the 
lattice. We find that the bulk gap is reduced firstly and 
the edge states are stable against a weak U, as shown in 
Figs. |4jc) to (f). When U is increased further, the bulk 
gap closes and the edge states disappear simultaneously. 
After that, a bulk gap with the character of the Mott gap 
occurs and no edge states reemerge anymore, as shown 
in Figs, ptg) and (h). Thus, we determine the phase 
boundary where the TBI disappears by a criteria that 
the bulk gap closes and the edge states disappear. Com- 
bining with the results described above, we can conclude 
that the TBI phase will make transition to the SL phase 
when A < 0.125i and to the easy-xy plane AF phase for 
A > 0.125i, as presented in the phase diagram of Fig. [2j 

According to the bulk-boundary correspondence pQ, 
the existence of gapless edge states depends on the topo- 
logical class of the bulk band structure. The transition 
from TBI (topologically nontrivial state) to MI (toplogi- 
cally trivial state) must undergo a gap closing process in 
the bulk. As far as we know, this process is demonstrated 
clearly for the first time by a systematic numerical cal- 
culation presented here. 

A first attempt to understand the evolution of the 
spectrums in the KMH model is to include the AF or- 
der parameter m A = —m B = \(n^ — (A and B 
denote the sublattice in Fig. [lja)) [17] . which is con- 
sidered as a result of the electron correlations. This 
gives rise to the mean field dispersion given by E(k) — 
±-y/e 2 (k) + (A — [/to/2) 2 , with e(k) the bare dispersion. 
When {/to/2 = A, the spin-orbit gap closes. Then an- 
other gap Um/2 — A with a character of the Mott gap 
opens up with the further increase of U. However, 
the evolution shown in Fig. [4] exhibits a more com- 
plex behavior, namely both the valence and conduction 
bands around K are split into two subbands. It im- 
plies that another interaction term is needed to be in- 
cluded. Comparing the numerical results for different U 
and A, we note that the band splitting around K de- 
pends on X 2 /U. This is the exchange integral in H 2 
coming from the second-order process of the spin-orbit 
interaction as described above. So, we rewrite Hi as [17] 



H \ = -(V 2 ) £««>) («W ~ a n a u)( a jt a ^ ~ aji.au) 
and choose another parameter x = ( a tf a jt ~ a li a jl) ■ By 
using m A and x as adjustable parameters, we can give a 
fit to the numerical results, which is plot as white dashed 
lines in Fig. [4] This simple fit provides a possible un- 
derstanding of the gap closing and reopening processes 
in the bulk. 

Finally, let us discuss the possible effect of electron 
correlations on the edge states. As shown in the inset 
of Fig. |4j we notice a visible reduction of the velocity in 
helical Dirac fermions at the edge in the TBI phase. This 
renormalization arising from the two-particle scattering 
between the left and right moving modes due to elec- 
tron correlations, which is allowed by the time reversal 
symmetry [27 1 135 ] . 

In summary, we have investigated the interplay be- 
tween the Hubbard interaction and the spin-orbit cou- 
pling in the Kane-Mele-Hubbard model with the varia- 
tional cluster approach. We map a detail U — A phase di- 
agram, in which the topological band insulator, the spin 
liquid, and the antiferromagnetic insulator arc identified. 
We have shown a nontrivial evolution of the bulk bands 
in the topological quantum phase transition. 

This work was supported by NSF-China and the 
MOST-China. XCX is also supported by US-DOE 
through the grant DE-FG02-04ER46124. 



[1] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 
(2010). 

[2] X. L. Qi and S. C. Zhang, |arXiv:1008.2026| (2010) . 

[3] C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 95, 146802 

(2005); C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 95, 

226801 (2005). 

[4] B. A. Bernevig, T. A. Hughes, and S. C. Zhang, Science 

314, 1757 (2006). 
[5] M. Konig, . Wiedmann, C. Briine, A. Roth, H. Buhmann, 

L. W. Molenkamp, X. L. Qi, and S. C. Zhang, Science 

318, 766 (2007). 
[6] D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. S. Hor, R. J. Cava, 

and M. Z. Hasan, Nature (London) 452, 970 (2008). 
[7] D. Hsieh, Y. Xia, D. Qian, L. Wray, J. H. Dil, F. Meier, 

J. Osterwalder, L. Patthey, J. G. Checkelsky, N. P. Ong, 

A. V. Fedorov, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, 

R. J. Cava, and M. Z. Hasan, Nature (London) 460, 1101 

(2009). 

[8] Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. 

Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. 

Hasan, Nat. Phys. 5, 398 (2009). 
[9] H. Zhang, C. X. Liu, X. L. Qi, X. Dai, Z. Fang, and S. 

-C. Zhang, Nat. Phys. 5, 438 (2009). 
[10] Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K. 

Mo, X. L. Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, 

S. C. Zhang, I. R. Fisher, Z. Hussain, and Z.-X. Shen, 

Science 325, 178 (2009). 
[11] D. Pesin and L. Balents, Nat. Phys. 6, 376 (2010). 
[12] A. Shitade, H. Katsura, J. Kune, X. -L Qi, S. -C. Zhang, 



5 



and N. Nagaosa, Phys. Rev. Lett. 102, 256403 (2009). 
[13] M. Liu, C. Chang, Z. Zhang, Y. Zhang, W. Ruan, K. He, 

L. -1. Wang, X. Chen, J. - F. Jia, S. -C. Zh ang, Q. -K. 

Xue, X. Ma, and Y. Wang, |arXTv:1011.1055| (2010). 
[14] S. Raghu, X. L. Qi, C. Honerkamp, and S. C. Zhang, 

Phys. Rev. Lett. 100, 156401 (2008). 
[15] K. Sun, H. Yao, E. Fradkin, and S. A. Kivelson, Phys. 

Rev. Lett. 103, 046811 (2009). 
[16] M. Dzero, K. Sun, V. Galitski, and P. Coleman, Phys. 

Rev. Lett. 104, 106408 (2010). 
[17] S. Rachel and K. L. Hur, Phys. Rev. B 82, 075106 (2010). 
[18] C. N. Varney, K. Sun, M. Rigol, and V. Galitski, Phys. 

Rev. B 82, 115125 (2010). 
[19] Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad, and A. 

Muramatsu, Nature (London) 464, 847 (2010). 
[20] M. Hohenadler T. C. Lang and F. F. Assaad, 

|arXiv:1011.5063|(2010). 



[21] D. Zheng, C. Wu, and G. M. Zhang, |arXiv:1011.5858 
(2010). 

[22] M. Potthoff, Eur. Phys. J. B 32, 429 (2003); M. Pot- 
thoff, M. Aichhorn, and C. Dahnken, Phys. Rev. Lett. 

91, 206402 (2003). 

[23] D. Senechal, P.-L. Lavertu, M.-A. Marois, and A.-M. S. 

Tremblay, Phys. Rev. Lett. 94, 156404 (2005). 
[24] M. Aichhorn, E. Arrigoni, M. Potthoff, and W. Hanke, 

Phys. Rev. B 74, 024508 (2006). 
[25] M. Balzer and M. Potthoff, Phys. Rev. B 82, 174441 

(2010). 

[26] D. Senechal, and A.-M. S. Tremblay, Phys. Rev. Lett. 

92, 126401 (2004). 

[27] C. Wu, B. A. Bernevig, and S. C. Zhang, Phys. Rev. Lett. 

96, 106401 (2006). 
[28] C. Xu and J.E. Moore, Phys. Rev. B 73, 045322 (2006) 



